Computing singularities of perturbation series 
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Many properties of current ab initio approaches to the quantum many-body problem, both pertur- 
bational or otherwise, are related to the singularity structure of Rayleigh-Schrodinger perturbation 
theory. A numerical procedure is presented that in principle computes the complete set of singu- 
larities, including the dominant singularity which limits the radius of convergence. The method 
approximates the singularities as eigenvalues of a certain generalized eigenvalue equation which is 
solved using iterative techniques. It relies on computation of the action of the perturbed Hamil- 
tonian on a vector, and does not rely on the terms in the perturbation series. Some illustrative 
model problems are studied, including a Helium-like model with 5-function interactions for which 
M0ller-Plesset perturbation theory is considered and the radius of convergence found. 
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I. INTRODUCTION 

Many-body perturbation theory (MBPT) has been 
one of the most popular approaches for ab initio many- 
body structure calculations, both in atomic, nuclear and 
chemical physics. Low-order M0ller-Plesset (MP) par- 
tial sums were for many years considered highly accurate 
and the method of choice for calculations of ground state 
energies. However, in recent years it has become clear 
that the convergence properties are not that simple, and 
that plain MBPT more often than not is divergent ■ 

Divergent series in this context can still be very useful. 
The series should be considered not as a final answer, 
but - as a Taylor series of a particular function with a 
particular singularity structure - be analyzed to obtain 
new ways of summing the series. Indeed the now ex- 
tremely popular coupled cluster method, which has to a 
large extent supplanted low-order MBPT as the most ef- 
fective method for ab initio structure calculations, can 
be described in terms of summations of selected classes 
of diagrams (i.e., selected terms in the series) to infinite 
order Q. Numerous other ways of resumming the series 
give improvements of the convergence, such as Pade or 
algebraic approximants @, [HHI • Especially in nuclear 
physics, summations of classes of diagrams to infinite or- 
der, such as the random-phase approximation jl2j . have 
wide-spread use. 

The performance of MBPT and the various resumma- 
tion techniques is determined by the singularity struc- 
ture of the energy eigenvalue maps E n (X), where A is 
the perturbation parameter. The determination of these 
singularities, which are of branch-point type, is therefore 
crucial, but also very involved. Current approaches use 
the terms in the series to estimate the location of the 
singularities, perhaps in combination with approximants 
[9l-tlll fl3Hla |. However, due to a theorem by Darboux 
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[ill [l4| the asymptotic form of the series only gives in- 
formation about the dominant singularity, i.e., the one 
closest to the origin, and such methods may also be sen- 
sitive to round-off errors [ll|. It is also possible to do 
(very expensive) parameter sweeps of A to locate avoided 
crossings |P6l - [l8| . thereby discovering empirically some 
singularities, but only those with small imaginary parts. 

In this article, we present a general and reliable nu- 
merical procedure for computing in principle the com- 
plete set of singularities of the eigenvalue maps E n (X) 
and a procedure for determining the dominant singular- 
ity in standard Rayleigh-Schrodinger (RS) perturbation 
theory from the results, thereby finding the radius of con- 
vergence (ROC) of the series. The method relies solely on 
being able to compute the action of the Hamiltonian on 
a vector, which is compatible with the common approach 
of using the full configuration- interaction (FCI) method- 
ology for computing the series terms [J, H, Ea, US] • We 
apply the numerical procedure to several examples and 
discuss the results. 

We have chosen the examples for their instructive na- 
ture and the fact that we can compare with an explicit 
analysis. We analyze a simple harmonic oscillator with 
a (5-function potential added plj j. a three-electron quan- 
tum wire model in one spatial dimension [22| , and an MP 
treatment of a Helium-like model with (5-function inter- 
actions, which was also considered recently in detail by 
Herman and Hagedorn [l8[ using parameter sweeps. In 
this paper we make conclusions about the ROC of this 
model. We use only very simple basis sets based on stan- 
dard discretization techniques. The two first examples il- 
lustrate the properties of our numerical procedure, while 
the final example illustrates an application of moderate 
complexity. 

Our method is based on the characterization of the 
sing ularities as branch points in the complex plane 
[l6l. l23l |24| . Those are equivalently the points A* where 
eigenvalues coalesce. It has been shown by the authors 
[251 ] that these points can be approximated to high pre- 
cision by solving a particular two-parameter eigenvalue 
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problem. More precisely, we find A(e) such that a pair 
eigenvalues have a small relative distance s, i.e., E n (X) 
and E m (X) = (1 + e)E n (X) are both eigenvalues. We 
adapt this result, exploit the structure of the Hamilto- 
nian matrix, and combine this with modern solvers for 
eigenvalue problems. 

The dominant singularity for RS perturbation theory 
for the ground state is the branch point A* o closest to 
the origin where E meets E n , n ^ [23|,[24|. The second 
part of the numerical method is a procedure that tracks 
the eigenvalue branches from the branch points A* to the 
origin, thereby determining if it is the dominant branch 
point A* ; o- We have chosen to focus on locating the dom- 
inant branch point since the ROC is a fundamental prop- 
erty of the perturbation series. The tracking procedure 
can equally well be applied to study other branch points. 

After discussing the numerical method in Section HH 
we apply it to the model problems in Section iHll Finally 
we present our conclusions in Section HVl 

II. METHOD 

A. Properties of RS perturbation series 

Consider a Hamiltonian matrix H of dimension N on 
the form 

H(X) = H a + XV, 

where V is treated as a perturbation, and where A is a 
complex parameter introduced for convenience. For the 
actual physical system we have A = A p h ys € R- Both Hq 
and V are Hcrmitian matrices. The eigenvalues E n (X) 
of H(X) are the N roots of the characteristic polynomial 
dct[H(X)—EI], where I is the identity matrix. The eigen- 
values E n (X) are the branches of an N- valued algebraic 
function, whose only singular points (denoted A* ) are in 
fact of branch-point type [H, [H, [24[ . 

For Hcrmitian matrices the branch-points come in 
complex conjugate pairs. There are no real branch 
points, and in the generic case (see Section III Dp all 
branch points are of square-root type. For sufficiently 
small A— A* the eigenvalues can be expanded in a Puiseux 
series around each branch point (24|. This is contained 
in Katz' theorem [26| which we state here: 

Theorem 1 Suppose H(X) = Hq + XV is generic in the 
sense that Hq and V are Hermitian and chosen at ran- 
dom. Then for any pair of branches E n and E m there ex- 
ists a branch point A* at which E n {X<,) = E m (X*) — b nm . 
Moreover, for sufficiently small A — A* there exists a con- 
stant c nm such that 

E n (X) = b nm + c nm (X - A,) 1 / 2 + 0(X - A,) 

and 

E m (X) = b nm - c nm (X - X,) 1 ' 2 + 0(X - A,), 



where it is to be understood that the same branch of the 
square-root function is to be used in both equations. 

Katz' theorem may be viewed as a generalization of the 
well-known Wigner-von Neumann non-crossing rule [24| . 
It is interesting that all eigenvalue pairs are involved 
at some branch point, which implies that the function 
Eo(X) actually can be analytically continued to any ex- 
cited state E n (X). 

Finite-dimensional Hamiltonians usually arise due to 
some discretization in form of a finite basis set, e.g., using 
the FCI methodology. The singularity structure of the 
full problem is richer than in the finite-dimensional case, 
but we postpone a brief discussion to Section IHI Al 

In RS perturbation theory for the ground state one 
computes a truncated Taylor series for Eg(X), viz, 

K 

E (X)^Y, E o,^ k + 0(X K+1 ), 

k=0 

which is an asymptotic series approximating Eo(X) as 
A — >• 0. The coefficients E$^ can be generated recur- 
sively by insertion into the eigenvalue problem for H(X) 
which gives a series usually represented in form of Feyn- 
man diagrams. The actual computation of the terms 
become increasingly complicated for higher-order terms 
for many-body systems (in practice, one rarely computes 
more than sixth-order series using diagrammatic tech- 
niques), but if -ff(A) is is available as a matrix or as 
a procedure that computes matrix-vector products, the 
high-order terms are straightforward to compute [1, [l9[ . 

One of the important questions we consider in this 
paper is whether the truncated series is convergent as 
K — > oo for A = Aphys, that is to say whether the ROC is 
greater than A p h ys or not. As a Taylor series, the ROC is 
given by (A^ol, where A*^ is the smallest branch point, 
called the dominant branch point, where the branch be- 
longing to £b(0) meets a different branch E n) n^ 0. We 
say that the E and E n branch at A*, [H, 113, S3, H • 
We remark that in other perturbation theories, like the 
folded diagram series for the effective interaction in nu- 
clear physics [2~ij , other branch points may be dominant. 

Thus, to compute the dominant singularity we are 
looking for the points A, £ C not on the real line such 
that E n (X*) = Eo(X^) for n ^ 0. The ROC is then 
R = nrin{|A»|} = |A* o|- Instrumental to this we consider 
the values of A such that the matrix H(X) = Ho + XV 
has a double eigenvalue. In what follows we call such a 
value of A a critical value. It usually is a branch point, 
but in some cases we also get spurious solutions. 

B. Computing the m smallest critical values 

We saw above that it is possible to characterize the sin- 
gularities of a perturbation series by computing A such 
that H(X) has a double eigenvalue. The problem of find- 
ing all A such that a matrix depending linearly on A has 
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a double eigenvalue has been considered by the authors 
elsewhere [25[ ■ The derivation of the method presented 
here is based on a result in Ref. [HJ stating that all solu- 
tions can be approximated by the solutions of a general- 
ized eigenvalue problem denned as follows. We let e > 
be a small scalar, called the regularization parameter, 
and define the matrices Ao(e), Ai(e) E C N xN by 

A (e) = -I®V + (l + s)V ®I 

and 

Ai(e) = I (g) H - (1 + e)H ® I, 

where <E> denotes the Kronecker product. Consider now 
the generalized eigenvalue problem 

AAo(e)w = Ai(e)v, (1) 

where v = v\ ® v-i. An important result is is that the 
solutions of (JTJ approximate all A such that H (A) has a 
double eigenvalue. 

Approximations of the eigenvalues E n (X) can be ob- 
tained from the eigenvalue problem for the matrix Hq + 
XV, where A satisfies ([I}. To shed a light on this ap- 
proach, one can prove that the eigenvalues of (QJ cor- 
respond to the set of values of A for which the matrix 
Hq + XV has two eigenvalues within a relative distance of 
e, i.e., a pair of eigenvalues of the form E and E(l + e). 
This is clearly a relaxation of the problem. A suitable 
choice of e, as well as other implementation aspects, have 
been studied in detail [25j|. This includes the exclusion 
of spurious solutions. It is also showed that the error in 
A behaves like 0(e 2 ). 

Although the above method allows us to compute all 
critical values of A, it may be computationally prohibitive 
for large problems. The computational complexity is de- 
termined by the solution of the generalized eigenvalue 
problem ([1]), which requires 0(N 6 ) operations if all eigen- 
values are computed with a general purpose method. To 
overcome this problem we use an iterative method known 
as the Arnoldi method [29j, to compute the m smallest 
eigenvalues of ([T]), where m is a given integer. 

The Arnoldi method generalizes the familiar Lanczos 
iterations employed in FCI calculations to non-Hermitian 
matrices, and only requires an efficient computation of 
the matrix-vector product associated with the eigenvalue 
problem. The matrix- vector product associated with ([T]) 
is 

y = A 1 (e)- 1 A (s)x. (2) 

Let X,Y G C NxN be such that x = vcc(X) and y = 

vec(F), where vec : C NxN — s- denotes the vectoriza- 
tion operation, i.e., stacking the columns of the matrix 
on top of each other. A key to the success of our method 
is that we can express the matrix-vector product © in 
terms of the matrices X and Y. By straightforward ma- 
nipulations using the rules of the Kronecker product we 
obtain 

H Q Y - (1 + e)YH T = -VX + (I + e)XV. (3) 



This matrix equation, where Y is the unknown, is a ma- 
trix equation known as a Sylvester equation. The right- 
hand side can be evaluated in 0(N 3 ) operations. The 
Sylvester equation can be solved in Q(N 3 ) operations by 
using the Bartels-Stewart algorithm [301 ] , which is a stan- 
dard method for Sylvester equations. Hence, by exploit- 
ing the structure in this way, the matrix-vector products 
of (JTJ can be efficiently computed in 0(N 3 ) operations. 
If Ho is diagonal, this can be improved to 0(N 2 ) opera- 
tions, which is seen as follows. 

Let Y =: (y\, . . . , j/at) and (c\, . . . , cn) the columns of 
the right-hand side of ©. Suppose the diagonal entries 
of Hq are i?o,i,ij i = 1, ■ ■ ■ , N . It is straighforward to 
show that column i of ([3]) can be written as the solution 
of a linear system with a diagonal matrix, 

yi=diag(di,...,d N )a (4) 

where 

dj = > , j = 1, • • • , N. 

"o,j,j ^ (1 + £)Ha.i,i 

Now note that we can compute a column vector of Y us- 
ing (Ql with only O(N) operations. Hence, the Sylvester 
equation corresponding to diagonal H can be solved in 
0{N 2 ) operations. In the simulations in Section MI CI 
we will use this a ppr oach, whereas using the Bartels- 
Stewart algorithm [30( turned out to be more robust in 
Section HLB 

Two additional properties of the Arnoldi method 
makes it particularly suitable for our purposes. 

• As we shall illustrate in Section Mil the overall al- 
gorithm for computing the dominant branch point, 
outlined in Section III CI typically requires only a 
small number of critical values m < JV. The 
Arnoldi method can be useful if not all eigenval- 
ues need to be computed. 

• If the chosen m is deemed insufficient, we wish to 
continue the iteration. The Arnoldi method can be 
easily resumed if more eigenvalues are needed. 

The procedure above describes a method which can be 
used for quite large systems since the complexity of the 
matrix- vector product is only 0(N 2 ). The matrices V 
and Hq stem from discretizations and we wish to be able 
to solve as large systems as possible. We will now use that 
for large problems (fine discretization/large basis set), a 
somewhat accurate guess is available by solving a corre- 
sponding smaller problem (coarser discretization/small 
basis set). 

Inverse iteration [29j is a method to compute one eigen- 
pair where a reasonable approximation of the eigenvalue 
is already available. Inverse iteration is, similar to the 
Arnoldi method, also only based on matrix vector prod- 
ucts. It however, does not involve any orthogonalization 
step. Since it is only based on matrix-vector products we 
can use ([3]) directly with inverse iteration. 



4 



By using the Arnoldi method for a coarse discretiza- 
tion, and inverse iteration for a finer discretization, we 
can solve very large problems in a reliable way in a multi- 
level fashion. 



C. Computing the dominant branch point 

The value A* ; o <G C is the hrst branch point of Eo(X). 
Since if(A».o) has a double eigenvalue, we can compute 
candidates for A*.o, i.e. the critical values, with the pro- 
cedure described in Section Hi B I It now remains to de- 
termine which one of the candidate solutions computed 
with the method in Section Hi Bl corresponds to A*,o- 

We will use a computational approach based on follow- 
ing paths from the candidate solution A to the origin. It 
is justified by the following technical result. 

Proposition 2 Consider a critical value A such that 
|A| < |A*,o|- Let p : [0, 1] — > C be a parametrization of a 
curve from p(0) = A to p(l) = such that \p(9)\ < |A| for 
9 £ [0, 1] . Assume that p does not pass directly through 
another critical value. Then two continuous eigenvalue 
functions [0, 1] 3 9 4 E n {6) and [0, 1] 3 4 E b (6) 
satisfying E a (9),E m (6) G a(H + p(9)V) for 9 <= [0, 1] 
and E n (0) = E m (0) are uniquely defined. Moreover, we 
have 

E n (l) ^ E o (0) and E m (l) + E o (0). 

Proof. The first statement follows from Rouche's The- 
orem [3l| . The second statement can be proven by con- 
tradiction. More precisely, the statement E n (l) = Eq(0) 
or E m (l) = Eq(0) contradicts with the assumption |A| < 
|A*,o|. □ 
From Proposition [2] it follows that A* : o is the smallest 
branch point for which one of the corresponding curves, 
E n or E m , terminates at Eq(0). Only the branch points 
with positive (or negative) imaginary parts are relevant 
and some may be spurious. Denoting all relevant numer- 
ical branch points by {Xk}^ =1 , where 

|Ai| < |A 2 | < ... < \X N ,\, 

this brings us to the following algorithm. 

Algorithm 3 (Computation of the ROC) 

1. Compute Eq(0), set i = 1. 

2. Consider Xi and continue the two corresponding 
branches E n and E m for 9 € [0, 1], i.e. from X = \i 
to X = 0. 

3. If one of the branches terminates at Eq(0), then 
stop 

else i = i + 1, go to step 2. 



We conclude this section with some implementation 
aspects. For step 2. critical values of the parameter A are 
needed in increasing magnitude. These can be computed 
by the Arnoldi algorithm as described in Section IIIBI 
The value of m is fixed before the iteration starts. If it 
turns out to be insufficient, the Arnoldi process can still 
be resumed as also outlined in Section IIIBI 

For the continuation process in step 3. we assume that 
curve p is linear, i.e. the curve corresponding to A^ satis- 
fies 

p(9) = (l-6)Xi. 

As the critical values are isolated points, this line does not 
contain other critical values, with probability one. For 
the continuation of the eigenvalues we follow the eigen- 
values by sampling the line between 9 G [0, 1] with suf- 
ficiently many points. In our applications, 21 sampling 
points were sufficient to follow the eigenvalues accurately 
and not dominate the computation time. 

Although we have chosen to to so in our implemen- 
tation, it is not necessary to compute the whole set of 
eigenvalues along p{9). Standard continuation techniques 
may instead be used, where continuity of the eigenvalues 
with respect to 9 is exploited. See for example the book 
by Seydel [3^ . 



D. A comment on genericity 
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As stated in Section lll Al the statements concerning the 
nature and location of the branch points A* of the eigen- 
value maps depend on the fact that H(X) is generic: a 
statement indicating that the matrix is a "typical" Her- 
mitian matrix. Statements about generic matrices hold 
with probability 1 when the matrix is chosen at random. 

On the other hand, Hamiltonians are rarely generic: 
symmetries such as angular momentum conservation or 
parity invariance lead to a natural block structure in 
H(X) so that the eigenvalue problem decouples into 
smaller, unrelated problems. It is easy to see that the 
non-crossing rule may be violated, and consequently that 
not all critical points of H(X) are branch points if not all 
symmetries are removed from the system. The numerical 
procedure then yields spurious real or complex solutions 
corresponding to violations of the non-crossing rule or 
the square-root branch point classification, respectively. 



III. NUMERICAL RESULTS 

A. Branch points in complete basis limit 

In this section we apply the numerical procedure to 
three model problems: a simple one-dimensional har- 
monic oscillator with a (5-function spike, which is exactly 
solvable |21[ and equivalent to a center-of-mass frame for- 
mulation of a parabolic two-electron quantum wire with 
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(5-function interactions, a three-electron parabolic quan- 
tum wire with smoothed Coulomb interactions [22j , and 
a helium-like model with (5-function nuclear and inter- 
electron interactions [l8|. In the latter example we con- 
sider M0ller-Plesset perturbation theory, while in the 
other examples we let the perturbation V be the bare 
inter-particle interactions. 

The Hamiltonian matrix is like in most MBPT ap- 
proaches an approximation of a partial differential op- 
erator H obtained by a finite basis expansion with dis- 
cretization parameter h, i.e., as h — > 0, the dimension 
N — > oo and the discrete spectrum approaches the ex- 
act limit under mild conditions. (Special care has to be 
taken for the continuum spectrum of H if it exists.) 

One may characterize the singularities of the eigen- 
value map of % as a or fj singularities (Toj . The a sin- 
gularities are complex-conjugate pairs of branch points 
with non-zero imaginary parts. These are also called 
"intruder states" [la ]. A finite-dimensional Hamiltonian 
only has a type branch points. The /? singularities are 
real branch points corresponding to a coalescence of an 
eigenvalue with the continuous spectrum. Baker argued 
this is a g eneric feature of unconfined fermion systems 
[3 HI] ■ The approximate Hamiltonian will, as long 
as it contains sufficiently good approximations to con- 
tinuum states, have a cluster of branch points near a /3 
singularity. As h — > 0, assuming that the basis set is 
in fact complete, the continuous spectrum is "filled out" 
with discrete points, and hence there will be many close 
crossings clustering around (but never equal to) a real 
value. 

To interpret and classify the numerically found A* we 
must consider the non-trivial limit h — > 0. In general, 
three cases can be expected: 

1. The branch point approaches a finite, complex 
value, and represents an a singularity. 

2. The branch point approaches infinity, in case of 
which a singularity of the perturbation series dis- 
appears for the exact Hamiltonian. 

3. The branch point approaches a finite real value. 
This can happen in two separate ways: (i) Since 
the branch points come in complex conjugate pairs 
this means that the limit actually becomes an an- 
alytic point as h — > 0, i.e., a violation of the non- 
crossing rule (which does not hold in the infinite- 
dimensional case), (ii) The real limit corresponds 
to a /3 singularity. In that case infinitely many 
branch points must approach the same real value 
as h -> (N -> oo). 



B. Harmonic oscillator with 5-function 

We consider the toy model Hamiltonian [2l| 
H(X) 



Any fixed A = A p h ys may be taken as the actual, physical 
value for the toy model. 

The eigenvalue problem (H — E)ip(x) = may be 
solved to arbitrary precision and represents a particularly 
instructive test-case for our numerical procedure. Parity 
symmetry allows us to focus on even eigenfunctions which 
includes the ground state. (The odd eigenfunctions are 
in fact trivial since S(x)ip(x) = in this case.) Figure Q] 
shows the eigenvalues of the even eigenfunctions as A is 
varied. 



2, 




FIG. 1: Eigenvalues of the harmonic oscillator with a 5- 
function. Only eigenvalues of even eigenfunctions are shown. 
Notice the crossing with the real axis around A ~ —0.6758, 
which will give rise to a spurious solution in the numerical 
method. 

Introducing the even-numbered harmonic oscillator ba- 
sis functions u n (x) = 4>2n(x) we obtain (Ho) nm = (2n + 
l/2)S nm and V nm = 4>2n(0)4>2m(0), the latter being a 
rank 1 matrix. Here, 



<j> n {x) = (2"n!V^r 1/2 ff„(z)e-* 2/2 , 



(5) 



W + 2 a +XS{X) - 



with H n (x) bein g th e standard Hermite polynomials. It 
has been shown [25l | that the numerical procedure will 
give spurious solutions of large magnitude since V has 
rank one; in this case |A*| ~ 10 12 . Also, one false real 
value arises for A such that H{\) has a zero eigenvalue, 
see Figure [TJ Note that all spurious solutions are easily 
detected. 

Figure [5] shows the smallest computed branch points 
for various N with the dominating A +j o inset. The results 
for the various TV indicate that the qualitative distribu- 
tion of branch points does not change much with the basis 
size. For iV -> oo we then estimate that RS perturbation 
theory will converge for all |A| < 2. 

We remark that it is not easy to find the branch points 
by doing a parameter sweep. Figure [1] does not reveal 
clear avoided crossings involving any pairs of eigenvalues, 
which is explained by the large imaginary parts of the 
various A*. 
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FIG. 2: The smallest branch points for the harmonic oscillator 
with a <5-function, computed for various matrix sizes N. The 
dominating branch point A«.o is shown inset. 
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FIG. 3: Eigenvalue branch paths as A is gradually decreased 
from A» to zero for two branch points of the harmonic os- 
cillator with ^-potential. Three eigenvalues are shown. The 
upper panel shows the paths for the dominant branch point 
A„,o, while the lower panel shows the paths for the first non- 
dominant branch point. 



We conclude this subsection by plotting the paths the 
eigenvalues trace out when A is gradually decreased from 
A* to zero, i.e., we consider X(9) = (1 — #)A* and plot 
eigenvalues as function of 8. Figure [3] shows the result 
for A„ o and one other branch point. This illustrates the 
continuation process described in Section III CI 



C. Three-electron quantum wire 

The next numerical calculation is on a one-dimensional 
model of a three-electron parabolic quantum wire, called 
so due to the quasi-one-dimensional confinement [22| . 
The electrons interact via a regularized Coulomb poten- 
tial of the form 

1 

u(x\, X2) oc 



X-2\ 



1 - 



where in our calculations we have set a = 0.1. The Hamil- 
tonian is then of the form 



1 d 2 

2<9x? 



1 



Xj), 



where we have introduced the parameter A which, in the 
chosen units, measures the relative strengths of the in- 
teractions compared to the semiconductor bulk and the 
size of the trap. Again, any fixed A = A p h ys can be taken 
to be the actual value. 

Due to the harmonic confinement, the spectrum of 
'H(A) is discrete for all A. It can be shown using a theorem 
due to Kato [23| that for all complex A all the eigenvalues 
depend analytically on A, i.e., there are no singularities at 
all. This is basically due to the boundedness of u(xi, x^)- 
Thus the ROC is infinite in the exact problem, and any 
perturbation approach should converge. A discretization 
will, however, necessarily produce branch point singular- 
ities, which will approach infinity or real values as the 
discretization is made finer. 

We use a standard discretization based on Slater de- 
terminants constructed from spin-orbitals on the form 
cf) n (x)xa{s), where 4> n (x) are the harmonic oscillator 
functions and x±i/2 are the spinor basis functions. 
For a given M we use all possible determinants created 
from spin-orbitals with y^.j rii < M. Thus, we include 
all unperturbed three-body harmonic oscillator states of 
energy less than M + 3/2. We restrict our attention to 
the lowest possible total spin projection S z = h. 

This yields matrices Ho and V of dimension N = 
0(M 2 ) when we separate out the center-of-mass mo- 
tion which is a dynamical symmetry - the center-of-mass 
moves like a free particle in a harmonic oscillator. We 
only consider even-parity wave functions, which includes 
the ground state. These are the only symmetries of the 
Hamiltonian operator, resulting in matrices for which the 
generic statements hold. 

Having obtained these matrices, we compute the 
branch points A* for e = 10~ 4 and also deduce the ROC 
using our numerical procedure. It is worthwhile to men- 
tion that in this case, the tracking procedure reveals that 
the dominant singularity is a critical value far from being 
the smallest. 

It is instructive to study the behavior of the ROC as 
function of M, as shown in Figure|3J As expected it seems 
to approach infinity linearly with M as the discretization 
becomes finer. 
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FIG. 4: The radius of convergence as function of the number 
of oscillator shells M for the quantum wire model. A clear 
linear tendency towards R — oo is shown. 



FIG. 6: The smallest branch points at M = 10 for the quan- 
tum wire model (circles) and paths (1 — 6>)A* used for deter- 
mining which branches E n (X) meet at A,. See also Figure [3 
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FIG. 5: The real part of the dominant singularity ReA*,o as 
function of M for the quantum wire model. It changes sign 
around M — 7, so that A*,o changes from a back-door to a 
front-door singularity. 
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FIG. 7: The eigenvalue Eo(X*) (squares) corresponding to 
the branch points A* in Figure [(J] Also shown are the paths 
E n [{\— 6>)A»] for the branching eigenvalues, determining which 
branches E n (X) actually meet at A,. 



If one tries to characterize A*^ as an intruder state, 
it is revealed in Figure [5] that RcA*.o changes sign as 
M is increased. This means that the characterization as 
"front-door" or "back-door intruders" [l6[ is dependent 
on the basis used. 

To illustrate the continuation procedure for determin- 
ing which eigenvalues branch a given A*, we have shown 
a number of the branch points involving the ground state 
E (X) and some other E n (X) in Figure [51 We construct 
a path X(8) = (1 — #)A* (also shown) and compute the 
eigenvalues of H(X(6)) that branch at A*. As these are 
continued to A(l) = 0, they will be equal to unperturbed 
energies, and the branches are easy to determine. In 



Fig. [7] the eigenvalue paths Eo(X(9)) and E n (X(6)) are 
shown. It is clearly seen that the ground state and some 
excited state meet at A* . This also shows how the contin- 
uation procedure may be used to find the non-dominant 
singularities involving the ground state, which can be 
taken as input for approximant construction. 



D. A helium-like model 

The final example is a helium-like model in one spa- 
tial dimension, also considered by Herman and Hagcdorn 
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1 1 81 ] - It shares many of the qualitative features with the 
true helium atom in three spatial dimensions, and our 
goal is to determine the ROC for a M0llcr-Plcsset cal- 
culation of the ground state energy. The model has the 
Hamiltonian 

= n a + V, 

where V = S(x± — x 2 ) is the inter-electron interaction. 
The two electrons also interact with the nucleus of charge 
Z via (5-function potentials. We set Z = 1.38 for the 
calculations. 

For M0ller-Plesset perturbation theory we rewrite the 
physical Hamiltonian as 

•H = -H +W HF + (V-W HF ), (6) 

where the Hartree-Fock operator W HF is defined in the 
usual way [l6| . Since we are going to treat V — W HF as a 
perturbation, we introduce a parameter A, viz, 

H HF (A) = Ho + W HF + A(V - U HF ). 

for which "H HF (1) = H. We now wish to determine 
whether the MP series converges, i.e., the radius of con- 
vergence must be greater than R = 1 . 

For this, we need to determine the Hartree-Fock basis 
and energies. We do this using a linear finite element 
basis over the interval [—L,L] using n sub- intervals of 
length Ax = 2L/n to obtain the usual Roothan-Hall 
equations which are solved iteratively [T(ij |. In our calcu- 
lations we have L = 15 and n = 1000, giving Aa; = 0.03. 
The exact solution to the Hartree-Fock ground state can 
be obtained in this case [l8| which provides a good check 
on the accuracy of the implementation. 

The discretization parameters are Ax, L and the num- 
ber of single-particle functions M we use in the MP series. 
However, we will consider the spatial discretization pa- 
rameters to be fixed and sufficient for an "exact" treat- 
ment of the one-body Hamiltonians, and instead only 
focus on M. 

The ground state is a singlet state, for which the spatial 
wave function is symmetric with respect to interchange 
of x\ and x%, and has positive parity, i.e., 

ip(xi,x 2 ) = 1p{X2,X 1 ) = ip(-xi,-x 2 ). 

There are otherwise no dynamical symmetries in this 
problem. 

Figure [8] shows a parameter sweep of the eigenvalues of 
% HF (A). In this case, it is clear that there is a back-door 
intruder around Re A* « —1.1, so i? < 1.1 is obvious. It is 
highly likely that it corresponds to a f3 type singularity 
in the exact problem, as there are clearly several close 
crossings (verified in our computations) with continuum 
states (with positive energy at A = 0) near A*. This is 
also supported by Herman and Hagedorn [jjjj . 




FIG. 8: (Color version online) Plot of the eigenvalues E n {\) 
for the helium-like model, using M = 42 single-particle func- 
tions. Note crossing around A ~ —1.1 with small imaginary 
part. We therefore expect interesting crossings with Ei and 
some other E n to have significant imaginary part since no 
other avoided crossings involving the ground state is visible. 



However, it is not clear whether or not there are other 
branch points with, say, a small real part and imaginary 
part ImA* = 0.7, which would imply R < 1 and hence 
an ultimate divergence of the MP series. 

The position or existence (5 singularities in the dis- 
crete problem is highly dependent on the basis chosen 
0, H3 1 so our conclusions must be taken relative to our 
discretization, which does include "diffuse" basis func- 
tions approximating the continuous spectrum. 
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FIG. 9: Convergence radius of the helium-like model as func- 
tion of the number of single-particle functions M. A rapid 
convergence towards R > 1.13 is seen. 

We compute, like in the preceding numerical examples, 
the dominant branch point A^o for various M and study 
its behavior. For M < 13 we used the Arnoldi method, 
while for M > 13 the inverse iteration method. It turns 
out that, in fact, R > 1; the avoided crossing does in- 
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FIG. 10: Estimated error in the ROC for the helium-like 
model as function of the number of single-particle functions 
M. An exponential error is observed, i.e., \R(M) — R(oo)\ ~ 
exp(-/3M). 



deed come from the dominant branch point. In Figure [9] 
the ROC as function of M is shown. Clearly, it stabi- 
lizes around some value R> 1.12. In Figure [TU] an error 
estimate is plotted, based on the largest M = 42, and 
clearly the error behaves like exp(— (3M) where (3 > is 
a constant. 

In Figure El the paths E Q (\(0)) and E 1 (\(6)) arc 
shown, where X(9) = (1 — 0)A*,o- This corresponds 
to Figure [7] for the quantum wire model. The path is 
surprisingly complicated, showing that the eigenvalues 
may take complicated deviations from their initial un- 
perturbed values as A varies in a straight line. 
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FIG. 11: Eigenvalue path for M = 42 for the dominant branch 
point A„,o- As A takes values along the path (1 — (9)A», the 
branching eigenvalues start at E n (X,) = E n t(X*) and moves 
to E n (0) and E n /(0), so that which branches meet can be 
determined. Here, n = and n' = 1. 



IV. CONCLUSION 



We have described a numerical procedure to determine 
the singularities of the eigenvalues E n {\) of Hq + XV. Us- 
ing a continuation technique that tracks eigenvalues as 
function of A, the dominant singularity can be found. A 
simple generalization of this will enable the classification 
of other singularities with respect to which eigenvalues 
branch at A*. By continuing Step 2 and Step 3 in Al- 
gorithm [3] after A* ; o has been found, one can find the 
secondary dominating branch point and so on. 

The method has been successfully applied to instruc- 
tive examples, and in particular the convergence of 
a M0ller-Plcsset perturbation series for a helium-like 
model was established. 

The most important virtue of the method is that it 
searches the whole complex plane for singularities, and 
also in principle can find all these. This allows a much 
more detailed mapping of the singularity structure than 
the standard methods based on the asymptotic form of 
the terms in the series. 

Computing the singularities of E n (X) is much harder 
than computing only E n (X p h ys )- One cannot hope to 
be able to compute the whole set of singularities for a 
very large many-body system. However, obtaining in- 
sight into the distribution of singularities for "typical" 
quantum systems, such as the examples considered in 
this paper and others [l7| . makes the construction and 
analysis of general resummation schemes easier [l(| H(| . 
Considering that popular approaches to the many-body 
problem such as coupled cluster methods can be viewed 
in terms of perturbation series only serves to emphasize 
the importance of calculations of singularity structures. 
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